####### SCRIPT FOR: 
####### SATISFACTION WITH DEMOCRACY IN AN ELECTORAL CYCLE

rm(list=ls())
gc()

setwd("C:/Users/miros/Desktop/Research/Impact of Individual-Level Factors on Electoral Systems/satisfaction with democracy")

#Disable scientific notations
options(scipen=999)

# Individual data
library(readxl)
idata <- read_excel("ess_data_r18.xlsx", sheet = "ind_data")
library(tidyr)
clean_idata <- drop_na(idata, i_date)
rm(idata)
#View(clean_idata)

#Coding dummies for electoral cycle periods
#Making sure the dates are specified as dates
clean_idata$i_date <- as.Date(clean_idata$i_date)
clean_idata$last_ele_plusd <- as.Date(clean_idata$last_ele_plusd)
clean_idata$next_ele <- as.Date(clean_idata$next_ele)
#Calculating relative duration between the last election and interview
clean_idata$elevssur_relative <- NA
clean_idata$elevssur_relative <- as.numeric(clean_idata$i_date - clean_idata$last_ele_plusd)/as.numeric(clean_idata$next_ele - clean_idata$last_ele_plusd)
#Adding dummies
clean_idata$elcyc_year <- NA
clean_idata$elcyc_year[clean_idata$elevssur_relative < 1.01] <- 1
clean_idata$elcyc_year[clean_idata$elevssur_relative < 0.875] <- 0.75
clean_idata$elcyc_year[clean_idata$elevssur_relative < 0.625] <- 0.5
clean_idata$elcyc_year[clean_idata$elevssur_relative < 0.375] <- 0.25
clean_idata$elcyc_year[clean_idata$elevssur_relative < 0.125] <- 0
clean_idata$elcyc_year <- as.factor(clean_idata$elcyc_year)
table(clean_idata$elcyc_year)

#Coding dummies for electoral cycle based on days (for robustness-check)
#Calculating number of days between the last election and interview
clean_idata$elevssur_days <- NA
clean_idata$elevssur_days <- as.numeric(clean_idata$i_date - clean_idata$last_ele_plusd)
#Adding dummies
clean_idata$elcyc_ddays <- NA
clean_idata$elcyc_ddays[clean_idata$elevssur_days < 1834] <- 1825
clean_idata$elcyc_ddays[clean_idata$elevssur_days < 1461] <- 1460
clean_idata$elcyc_ddays[clean_idata$elevssur_days < 1278] <- 1095
clean_idata$elcyc_ddays[clean_idata$elevssur_days <  913] <-  730
clean_idata$elcyc_ddays[clean_idata$elevssur_days <  548] <-  365
clean_idata$elcyc_ddays[clean_idata$elevssur_days <  183] <-    0
clean_idata$elcyc_ddays <- as.factor(clean_idata$elcyc_ddays)
table(clean_idata$elcyc_ddays)

#Changing label for winners/losers
clean_idata$winner[clean_idata$winner == 1] <- "Winner"
clean_idata$winner[clean_idata$winner == 0] <- "Loser"
clean_idata$winner[clean_idata$winner == 99] <- NA
table(clean_idata$winner)

#Turning absenteers into a separate category
#(not used in the analysis)
#clean_idata$winner1 <- NA
#clean_idata$winner1[is.na(clean_idata$winner)] <- "Absenter"
#clean_idata$winner1[clean_idata$winner == 1] <- "Winner"
#clean_idata$winner1[clean_idata$winner == 0] <- "Loser"
#table(clean_idata$winner1)




#FIRST PART OF THE ANALYSIS: ELECTORAL SYSTEMS
#Dividing countries into electoral system categories
clean_idata$el_sys <- NA
clean_idata$el_sys[clean_idata$cntry == "FR" | clean_idata$cntry == "GB"]                             <- "Majoritarian system"
clean_idata$el_sys[clean_idata$cntry == "DE" | clean_idata$cntry == "GR" | clean_idata$cntry == "HU"| 
                   clean_idata$cntry == "IT" | clean_idata$cntry == "LT" | clean_idata$cntry == "RU"] <- "Mixed system"
clean_idata$el_sys[clean_idata$cntry == "UA" & clean_idata$el_year == 2002]                           <- "Mixed system"
clean_idata$el_sys[clean_idata$cntry == "UA" & clean_idata$el_year == 2012]                           <- "Mixed system"
clean_idata$el_sys[clean_idata$cntry == "AT" | clean_idata$cntry == "BE" | clean_idata$cntry == "BG"| 
                   clean_idata$cntry == "HR" | clean_idata$cntry == "CY" | clean_idata$cntry == "CZ"| 
                   clean_idata$cntry == "DK" | clean_idata$cntry == "EE" | clean_idata$cntry == "FI"| 
                   clean_idata$cntry == "IS" | clean_idata$cntry == "IE" | clean_idata$cntry == "IL"| 
                   clean_idata$cntry == "LU" | clean_idata$cntry == "NL" | clean_idata$cntry == "NO"| 
                   clean_idata$cntry == "PL" | clean_idata$cntry == "PT" | clean_idata$cntry == "SK"| 
                   clean_idata$cntry == "SI" | clean_idata$cntry == "ES" | clean_idata$cntry == "SE"| 
                   clean_idata$cntry == "CH" | clean_idata$cntry == "TR"]                             <- "Proportional system"
clean_idata$el_sys[clean_idata$cntry == "UA" & clean_idata$el_year == 2006]                           <- "Proportional system"
clean_idata$el_sys[clean_idata$cntry == "UA" & clean_idata$el_year == 2007]                           <- "Proportional system"

#Ukraine implemented several electoral reforms, thus it is split
clean_idata$cntry[clean_idata$cntry == "UA" & clean_idata$el_year == 2006] <- "UA (2006, 2007)"
clean_idata$cntry[clean_idata$cntry == "UA" & clean_idata$el_year == 2007] <- "UA (2006, 2007)"
clean_idata$cntry[clean_idata$cntry == "UA" & clean_idata$el_year == 2002] <- "UA (2002)"

clean_idata$country[clean_idata$country == "Ukraine" & clean_idata$el_year == 2006] <- "Ukraine (2006, 2007)"
clean_idata$country[clean_idata$country == "Ukraine" & clean_idata$el_year == 2007] <- "Ukraine (2006, 2007)"
clean_idata$country[clean_idata$country == "Ukraine" & clean_idata$el_year == 2002] <- "Ukraine (2002)"



#Creating new data frame for Models 1 to 5 to sync their number of observations
data_m15 <- data.frame(clean_idata$cntry, clean_idata$country, clean_idata$essround, clean_idata$p_ess_survey, clean_idata$el_cyc, clean_idata$el_sys,
                       clean_idata$stfdem, clean_idata$elcyc_year, clean_idata$elcyc_ddays, clean_idata$winner, clean_idata$LSq,
                       clean_idata$e_democ, clean_idata$length_e_democ, clean_idata$snap_el, clean_idata$polarization,
                       clean_idata$agea, clean_idata$edulvla, clean_idata$female, clean_idata$stfeco, clean_idata$contplt, 
                       clean_idata$pspwght)
data_m15 <- na.omit(data_m15)
colnames(data_m15) <- c("cntry", "country", "essround", "p_ess_survey", "el_cyc", "el_sys",
                        "stf_dem", "elcyc_year", "elcyc_ddays", "winner", "LSq",
                        "dem_dev", "dem_trad", "snap_el", "polarization", 
                        "age", "edu_narrow", "female", "stfeco", "contplt",
                        "pspwght")
#write.csv(data_m15, "Cases used for electoral system analysis.csv")
summary(data_m15)

#Calculating averages for electoral system categories to display them
library(stats)
a.cntry.avg   <- aggregate(data_m15$LSq, list(data_m15$country), mean)
b.cntry.elsys <- unique(data.frame(data_m15$el_sys, data_m15$country))
c.elsys.avg   <- aggregate(data_m15$LSq, list(data_m15$el_sys), mean)

elsys.avg <- merge(merge(a.cntry.avg, b.cntry.elsys, by.x = "Group.1", by.y = "data_m15.country"), c.elsys.avg, by.x = "data_m15.el_sys", by.y = "Group.1")
colnames(elsys.avg) <- c("Electoral System Category", "Country", "Average Country Disproportionality", "Average Electoral System Disproportionality")
elsys.avg <- elsys.avg[c(1, 4, 2, 3)]
#Table A2: Average electoral system disproportionality
#write.csv(elsys.avg, "Average proportionality of electoral systems.csv")
rm(a.cntry.avg, b.cntry.elsys, c.elsys.avg)


#Model 1: null-model
library(lme4)
library(AICcmodavg)
library(performance)
ind_SWD.1 <- lmer(stf_dem ~ 1 + (1 | cntry/p_ess_survey),
                  data=data_m15, REML=FALSE, weights=pspwght) 
summary(ind_SWD.1)
icc(ind_SWD.1)

#Model 2: Winners & electoral cycle
ind_SWD.2 <- lmer(stf_dem ~ elcyc_year*winner
                  + (1 | cntry/p_ess_survey),
                  data=data_m15, REML=FALSE, weights=pspwght) 
summary(ind_SWD.2)
icc(ind_SWD.2)

#Model 3: Full interaction, no controls
ind_SWD.3 <- lmer(stf_dem ~ elcyc_year*winner*LSq
                  + (1 | cntry/p_ess_survey),
                  data=data_m15, REML=FALSE, weights=pspwght) 
summary(ind_SWD.3)
icc(ind_SWD.3)

#Model 4: Winners + electoral cycle periods + controls
ind_SWD.4 <- lmer(stf_dem ~ elcyc_year*winner
                  + dem_dev + dem_trad + snap_el + polarization
                  + age + edu_narrow + female  
                  + (1 | cntry/p_ess_survey),
                  data=data_m15, REML=FALSE, weights=pspwght) 
summary(ind_SWD.4)
icc(ind_SWD.4)

#Model 5: Full model
ind_SWD.5 <- lmer(stf_dem ~ elcyc_year*winner*LSq
                  + dem_dev + dem_trad + snap_el + polarization
                  + age + edu_narrow + female
                  + (1 | cntry/p_ess_survey) ,
                  data=data_m15, REML=FALSE, weights=pspwght) 
summary(ind_SWD.5)
icc(ind_SWD.5)

#Table 1: Output with models 1-5
library(stargazer)
stargazer(ind_SWD.1, ind_SWD.2, ind_SWD.3, ind_SWD.4, ind_SWD.5, 
          type="html", out="model1_5.htm", 
          single.row = TRUE, star.cutoffs = c(0.05, 0.01, 0.001), digits = 3)


#Figure 5: Electoral systems
library(ggplot2)
library(sjPlot)
library(sjmisc)
#png(file = "Fig5_electoral_systems.png", width = 6150, height = 2650, res = 600)
ind_SWD.5p <- plot_model(ind_SWD.5,
                         type = "pred",
                         terms = c("elcyc_year", "winner", "LSq [4.11, 7.35, 17.01]"),
                         title = "",
                         colors = "gs") + 
  #geom_line(stat="smooth", method="loess") +
  coord_cartesian(xlim = c(0.9, 5.1), ylim = c(3, 7), expand = FALSE) +
  scale_y_continuous(breaks = c(3, 3.5, 4, 4.5, 5, 5.5, 6, 6.5, 7),
                     labels = c("3", "", "4", "", "5", "", "6", "", "7")) +
  labs(x = "Proportion of electoral cycle time elapsed",
       y = "Satisfaction with democracy \n[0 = Extremely dissatisfied; 1 = Extremely satisfied]") +
  theme(plot.title = element_blank(),
        legend.position="bottom",
        legend.title=element_blank(),
        legend.background = element_rect(colour = 'black', fill = 'white', linetype='solid'),
        legend.key=element_blank(),
        panel.border = element_rect(fill = "transparent"),
        strip.background = element_rect(colour = "black"),
        panel.grid.major.x = element_line(colour = "gray95"),
        panel.grid.minor.x = element_blank(),
        panel.grid.major.y = element_line(colour = "gray95"),
        panel.grid.minor.y = element_line(colour = "gray95"),
        panel.background = element_blank())
#changing the facet labels
ind_SWD.5p$data$facet <- factor(ind_SWD.5p$data$facet, 
                                levels = c("LSq = 4.11", "LSq = 7.35", "LSq = 17.01"),
                                labels = c("Proportional systems (disproportionality = 4.11)", 
                                           "Mixed systems (disproportionality = 7.35)", 
                                           "Majoritarian systems (disproportionality = 17.01)"))
ind_SWD.5p
#dev.off()



#ROBUSTNESS-CHECK for the first part of the analysis
#Days instead of periods on electoral cycle
library(lme4)
library(AICcmodavg)
library(performance)
#Model A1 (Robustness-check): null-model
ind_SWD.A1 <- lmer(stf_dem ~ 1 + (1 | cntry/p_ess_survey),
                  data=data_m15, REML=FALSE, weights=pspwght) 
summary(ind_SWD.A1)
icc(ind_SWD.A1)

#Model A2 (Robustness-check): Winners & electoral cycle
ind_SWD.A2 <- lmer(stf_dem ~ elcyc_ddays*winner
                  + (1 | cntry/p_ess_survey),
                  data=data_m15, REML=FALSE, weights=pspwght) 
summary(ind_SWD.A2)
icc(ind_SWD.A2)

#Model A3 (Robustness-check): Full interaction, no controls
ind_SWD.A3 <- lmer(stf_dem ~ elcyc_ddays*winner*LSq
                  + (1 | cntry/p_ess_survey),
                  data=data_m15, REML=FALSE, weights=pspwght) 
summary(ind_SWD.A3)
icc(ind_SWD.A3)

#Model A4 (Robustness-check): Winners + electoral cycle periods + controls
ind_SWD.A4 <- lmer(stf_dem ~ elcyc_ddays*winner
                  + dem_dev + dem_trad + snap_el + polarization
                  + age + edu_narrow + female
                  + (1 | cntry/p_ess_survey),
                  data=data_m15, REML=FALSE, weights=pspwght) 
summary(ind_SWD.A4)
icc(ind_SWD.A4)

#Model A5 (Robustness-check): Full model
ind_SWD.A5 <- lmer(stf_dem ~ elcyc_ddays*winner*LSq
                  + dem_dev + dem_trad + snap_el + polarization
                  + age + edu_narrow + female 
                  + (1 | cntry/p_ess_survey) ,
                  data=data_m15, REML=FALSE, weights=pspwght) 
summary(ind_SWD.A5)
icc(ind_SWD.A5)

#Table A4: Output with models A1-A5
library(stargazer)
stargazer(ind_SWD.A1, ind_SWD.A2, ind_SWD.A3, ind_SWD.A4, ind_SWD.A5, 
          type="html", out="modelA1_A5_robustness_check.htm", 
          single.row = TRUE, star.cutoffs = c(0.05, 0.01, 0.001), digits = 3)

#Figure A1: Electoral systems
library(ggplot2)
library(sjPlot)
library(sjmisc)
#png(file = "FigA1_electoral_systems.png", width = 6150, height = 2650, res = 600)
ind_SWD.A5p <- plot_model(ind_SWD.A5,
                         type = "pred",
                         terms = c("elcyc_ddays", "winner", "LSq [4.11, 7.35, 17.01]"),
                         title = "",
                         colors = "gs") + 
  #geom_line(stat="smooth", method="loess") +
  coord_cartesian(xlim = c(0.9, 6.1), ylim = c(3, 7), expand = FALSE) +
  scale_x_continuous(breaks = c(1:6), labels = c(0:5)) +
  scale_y_continuous(breaks = c(3, 3.5, 4, 4.5, 5, 5.5, 6, 6.5, 7),
                     labels = c("3", "", "4", "", "5", "", "6", "", "7")) +
  labs(x = "Proportion of electoral cycle time elapsed",
       y = "Satisfaction with democracy \n[0 = Extremely dissatisfied; 1 = Extremely satisfied]") +
  theme(plot.title = element_blank(),
        legend.position="bottom",
        legend.title=element_blank(),
        legend.background = element_rect(colour = 'black', fill = 'white', linetype='solid'),
        legend.key=element_blank(),
        panel.border = element_rect(fill = "transparent"),
        strip.background = element_rect(colour = "black"),
        panel.grid.major.x = element_line(colour = "gray95"),
        panel.grid.minor.x = element_blank(),
        panel.grid.major.y = element_line(colour = "gray95"),
        panel.grid.minor.y = element_line(colour = "gray95"),
        panel.background = element_blank())
#changing the facet labels
ind_SWD.A5p$data$facet <- factor(ind_SWD.A5p$data$facet, 
                                levels = c("LSq = 4.11", "LSq = 7.35", "LSq = 17.01"),
                                labels = c("Proportional systems (disproportionality = 4.11)", 
                                           "Mixed systems (disproportionality = 7.35)", 
                                           "Majoritarian systems (disproportionality = 17.01)"))
ind_SWD.A5p
#dev.off()





#SECOND PART OF THE ANALYSIS: PREDICTABILITY OF PARTY SYSTEMS
#Putting Ukraine back together
clean_idata$cntry[clean_idata$cntry == "UA (2006, 2007)" | clean_idata$cntry == "UA (2002)"] <- "UA"
clean_idata$country[clean_idata$country == "Ukraine (2006, 2007)" | clean_idata$country == "Ukraine (2002)"] <- "Ukraine"

#Creating a new data frame for Models 6 to 10 to sync their number of observations
data_m610 <- data.frame(clean_idata$cntry, clean_idata$country, clean_idata$essround, clean_idata$p_ess_survey, clean_idata$el_cyc, 
                       clean_idata$stfdem, clean_idata$elcyc_year, clean_idata$elcyc_ddays, clean_idata$winner, clean_idata$closure,
                       clean_idata$e_democ, clean_idata$length_e_democ, clean_idata$snap_el, clean_idata$polarization,
                       clean_idata$agea, clean_idata$edulvla, clean_idata$female, clean_idata$stfeco, clean_idata$contplt, 
                       clean_idata$pspwght)
data_m610 <- na.omit(data_m610)
colnames(data_m610) <- c("cntry", "country", "essround", "p_ess_survey", "el_cyc", 
                        "stf_dem", "elcyc_year", "elcyc_ddays", "winner", "closure",
                        "dem_dev", "dem_trad", "snap_el", "polarization", "stfeco", "contplt", 
                        "age", "edu_narrow", "female", 
                        "pspwght")
summary(data_m610)
#write.csv(data_m610, "Cases used for predictability analysis.csv")


#Calculating country averages for party system closure
library(stats)
cntry.cls.avg   <- aggregate(data_m610$closure, list(data_m610$cntry, data_m610$country), mean)
library(dplyr)
cntry.cls.avg <- arrange(cntry.cls.avg, x)
colnames(cntry.cls.avg) <- c("Abbreviation", "Country", "Average Closure")
#Table A3: Average electoral system disproportionality
#write.csv(cntry.cls.avg, "Average closure.csv")

#Model 6: null-model
library(lme4)
library(AICcmodavg)
library(performance)
ind_SWD.6 <- lmer(stf_dem ~ 1 + (1 | cntry/p_ess_survey),
                  data=data_m610, REML=FALSE, weights=pspwght) 
summary(ind_SWD.6)
icc(ind_SWD.6)

#Model 7: Winners & electoral cycle
ind_SWD.7 <- lmer(stf_dem ~ elcyc_year*winner
                  + (1 | cntry/p_ess_survey),
                  data=data_m610, REML=FALSE, weights=pspwght) 
summary(ind_SWD.7)
icc(ind_SWD.7)

#Model 8: Full interaction, no controls
ind_SWD.8 <- lmer(stf_dem ~ elcyc_year*winner*closure
                  + (1 | cntry/p_ess_survey),
                  data=data_m610, REML=FALSE, weights=pspwght) 
summary(ind_SWD.8)
icc(ind_SWD.8)

#Model 9: Winners + electoral cycle periods + controls
ind_SWD.9 <- lmer(stf_dem ~ elcyc_year*winner
                  + dem_dev + dem_trad + snap_el + polarization
                  + age + edu_narrow + female
                  + (1 | cntry/p_ess_survey),
                  data=data_m610, REML=FALSE, weights=pspwght) 
summary(ind_SWD.9)
icc(ind_SWD.9)

#Model 10: Full model
ind_SWD.10 <- lmer(stf_dem ~ elcyc_year*winner*closure
                  + dem_dev + dem_trad + snap_el + polarization
                  + age + edu_narrow + female
                  + (1 | cntry/p_ess_survey) ,
                  data=data_m610, REML=FALSE, weights=pspwght) 
summary(ind_SWD.10)
icc(ind_SWD.10)

#Table 2: Output with models 6-10
library(stargazer)
stargazer(ind_SWD.6, ind_SWD.7, ind_SWD.8, ind_SWD.9, ind_SWD.10, 
          type="html", out="model6_10.htm", 
          single.row = TRUE, star.cutoffs = c(0.05, 0.01, 0.001))


#Figure 6: Predictability of the political development
library(ggplot2)
library(sjPlot)
library(sjmisc)
#png(file = "Fig6_predictability.png", width = 6150, height = 2650, res = 600)
ind_SWD.10p <- plot_model(ind_SWD.10,
                         type = "pred",
                         terms = c("elcyc_year", "winner", "closure [70.34, 90.70, 98.99]"),
                         title = "",
                         colors = "gs") +
  #geom_line(stat="smooth", method="loess") +
  coord_cartesian(xlim = c(0.9, 5.1), ylim = c(3, 7), expand = FALSE) +
  scale_y_continuous(breaks = c(3, 3.5, 4, 4.5, 5, 5.5, 6, 6.5, 7),
                     labels = c("3", "", "4", "", "5", "", "6", "", "7")) +
  labs(x = "Proportion of electoral cycle time elapsed",
       y = "Satisfaction with democracy \n[0 = Extremely dissatisfied; 1 = Extremely satisfied]") +
  theme(plot.title = element_blank(),
        legend.position="bottom",
        legend.title=element_blank(),
        legend.background = element_rect(colour = 'black', fill = 'white', linetype='solid'),
        legend.key=element_blank(),
        panel.border = element_rect(fill = "transparent"),
        strip.background = element_rect(colour = "black"),
        panel.grid.major.x = element_line(colour = "gray95"),
        panel.grid.minor.x = element_blank(),
        panel.grid.major.y = element_line(colour = "gray95"),
        panel.grid.minor.y = element_line(colour = "gray95"),
        panel.background = element_blank())
#changing and ordering the facet labels
ind_SWD.10p$data$facet <- factor(ind_SWD.10p$data$facet, 
                                 levels = c("closure = 70.34", "closure = 90.7", "closure = 98.99"),
                                 labels = c("Low predictability (closure = 70.34; avg for Ukraine)", 
                                            "Medium predictability (closure = 90.70; avg for Belgium)", 
                                            "High predictability (closure = 98.99; avg for Switzerland)"))
ind_SWD.10p
#dev.off()



#ROBUSTNESS-CHECK for the second part of the analysis
#Days instead of periods on electoral cycle
library(lme4)
library(AICcmodavg)
#Model A6: null-model
ind_SWD.A6 <- lmer(stf_dem ~ 1 + (1 | cntry/p_ess_survey),
                  data=data_m610, REML=FALSE, weights=pspwght) 
summary(ind_SWD.A6)
icc(ind_SWD.A6)

#Model A7: Winners & electoral cycle
ind_SWD.A7 <- lmer(stf_dem ~ elcyc_ddays*winner
                  + (1 | cntry/p_ess_survey),
                  data=data_m610, REML=FALSE, weights=pspwght) 
summary(ind_SWD.A7)
icc(ind_SWD.A7)

#Model A8: Full interaction, no controls
ind_SWD.A8 <- lmer(stf_dem ~ elcyc_ddays*winner*closure
                  + (1 | cntry/p_ess_survey),
                  data=data_m610, REML=FALSE, weights=pspwght) 
summary(ind_SWD.A8)
icc(ind_SWD.A8)

#Model A9: Winners + electoral cycle periods + controls
ind_SWD.A9 <- lmer(stf_dem ~ elcyc_ddays*winner
                  + dem_dev + dem_trad + snap_el + polarization
                  + age + edu_narrow + female 
                  + (1 | cntry/p_ess_survey),
                  data=data_m610, REML=FALSE, weights=pspwght) 
summary(ind_SWD.A9)
icc(ind_SWD.A9)

#Model 10: Full model
ind_SWD.A10 <- lmer(stf_dem ~ elcyc_ddays*winner*closure
                   + dem_dev + dem_trad + snap_el + polarization
                   + age + edu_narrow + female 
                   + (1 | cntry/p_ess_survey) ,
                   data=data_m610, REML=FALSE, weights=pspwght) 
summary(ind_SWD.A10)
icc(ind_SWD.A10)


#Table A5: Output for comparison of the models 4-6 with their respecified versions
library(stargazer)
stargazer(ind_SWD.A6, ind_SWD.A7, ind_SWD.A8, ind_SWD.A9, ind_SWD.A10, 
          type="html", out="modelA6_A10_robustness_check.htm", 
          single.row = TRUE, star.cutoffs = c(0.05, 0.01, 0.001))

#Figure A2: Predictability of the political development
library(ggplot2)
library(sjPlot)
library(sjmisc)
#png(file = "FigA2_predictability.png", width = 6150, height = 2650, res = 600)
ind_SWD.A10p <- plot_model(ind_SWD.A10,
                          type = "pred",
                          terms = c("elcyc_ddays", "winner", "closure [70.34, 90.70, 98.99]"),
                          title = "",
                          colors = "gs") +
  #geom_line(stat="smooth", method="loess") +
  coord_cartesian(xlim = c(0.9, 6.1), ylim = c(3, 7), expand = FALSE) +
  scale_x_continuous(breaks = c(1:6), labels = c(0:5)) +
  scale_y_continuous(breaks = c(3, 3.5, 4, 4.5, 5, 5.5, 6, 6.5, 7),
                     labels = c("3", "", "4", "", "5", "", "6", "", "7")) +
  labs(x = "Proportion of electoral cycle time elapsed",
       y = "Satisfaction with democracy \n[0 = Extremely dissatisfied; 1 = Extremely satisfied]") +
  theme(plot.title = element_blank(),
        legend.position="bottom",
        legend.title=element_blank(),
        legend.background = element_rect(colour = 'black', fill = 'white', linetype='solid'),
        legend.key=element_blank(),
        panel.border = element_rect(fill = "transparent"),
        strip.background = element_rect(colour = "black"),
        panel.grid.major.x = element_line(colour = "gray95"),
        panel.grid.minor.x = element_blank(),
        panel.grid.major.y = element_line(colour = "gray95"),
        panel.grid.minor.y = element_line(colour = "gray95"),
        panel.background = element_blank())
#changing and ordering the facet labels
ind_SWD.A10p$data$facet <- factor(ind_SWD.A10p$data$facet, 
                                 levels = c("closure = 70.34", "closure = 90.7", "closure = 98.99"),
                                 labels = c("Low predictability (closure = 70.34; avg for Ukraine)", 
                                            "Medium predictability (closure = 90.70; avg for Belgium)", 
                                            "High predictability (closure = 98.99; avg for Switzerland)"))
ind_SWD.A10p
#dev.off()